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Abstract. We study the photoionization and autoionization of Helium atom subject 
to ultrashort laser pulses by using a Feshbach formalism in the time domain. We 
solve the time-dependent Schrodinger equation in terms of a configuration interaction 
(CI) spectral method, in which the total wavefunction is expanded with configurations 
defined within bound-like (Q) and scattering- like (V) halfspaces. The method allows 
^ ■ one to provide accurate descriptions of both the atomic structure (energy positions 

and widths) and the photodynamics. We illustrate our approach by i) calculating the 
time-resolved one-photon ionization below the He + (n=2) ionization threshold, from 
l 1 S e and 2 1 P° initial states, then reaching the lowest autoionizing states of 1 S e , 1 P° 
and 1 D e final symmetries ii) studing the temporal formation of the Fano profile of 1 P° 
resonances and iii) showing its performance in obtaining the perturbative long-time 
limit of one- and two-photon ionization cross sections using ultrashort laser pulses 
following a recently developed procedure in Phys. Rev. A, 77, 032716 (2008). 
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1. Introduction 



In the last decade the experimental development of high intensity, high frequency 
ultrashort laser pulses down to the attosecond time scale has led to the investigation 
of electron photodynamics in its natural time scale [TJ [2]. In general, the 
temporal behaviour of many-electron atoms subject to a strong laser field still 
remains a challenging theoretical problem when dealing with multiple excitation and 
ionization, including Auger phenomena, which requires to account for the appropriate 
representation of the electron correlation (responsible for the autoionization) along with 
the solution of the time-dependent Schrodinger equation. For instance, it is interesting 
to study the coherent excitation and decay of metastable superexcited states using 
pump-probe laser pulses with time delays comparable in magnitude to autoionization 
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lifetimes in order to trace the fast dynamics involved. The decay dynamics of doubly 
excited states in the Helium atom has already been addressed using time- dependent 
close-coupling methods implemented in a numerical lattice [3]. Also, attosecond pump 
probe laser schemes have been proposed to probing ultrafast electron motion in singly 
excited jl] and doubly excited states in Helium [3J. In particular, the latter work is 
based on a sophisticated multichannel close-coupling method using B-splines with very 
large radial boxes, large L-values in the angular momentum expansion, introduction 
of absorbing potentials, etc., ingredients that require in general a high computational 
effort. 

Considering fast evolving phenomena, it has been of recent interest the time- 
resolved formation of Fano profiles in the atomic photoionization spectra, using both 
simplified models and ab initio methods [El [7J [8]. In this respect, experimental 
attosecond resolution of the He autoionization process has been recently reported 
[9], demonstrating control over the two interfering paths, direct ionization and 
autoionization, that shape the profile of the Fano resonance pU] (for a recent study 
on the modification of the asymmetry in the Fano profiles with ultrashort pulses see 
also [llj). New emergent phenomena arising in the autoionization dynamics of laser- 
dressed atoms (pump probe schemes) where excited resonances are coupled by an IR 
field, have been mostly studied using simplified models [121 02] , that although reproduce 
the relevant features of experiments and support the physical explanation, they may 
prove insufficient in other scenarios that require fully ab initio methods [14]. 

In this work we describe a time dependent Feshbach method applied to atomic 
photoionization of He subject to ultrashort laser pulses. The time independent Feshbach 
formalism has been widely applied in atomic physics in the last decades but the time 
dependent formalism in the present form has been used only in the molecular context 
recently [13]. A similar time- dependent approach (but not identical in the form of 
the wave packet expansion and the required final projections) was proposed to study 
the resonant and non resonant ionization of He by XUV intense ultrashort laser pulses 
[16] . In particular, we describe here some practical details to generate discretized non 
resonant continuum states for a given selected energy from the solution of the eigenvalue 
problem. To gauge the performance of the method, we provide simple illustrations like 
the computation of one- and two-photon ionization cross section in He. For weak fields, 
the dependence of the one-photon ionization probability with the pulse duration obeys 
only to spectral effects. In fact, the temporal dependence in the transition amplitudes 
calculated in time dependent first order perturbation theory is exactly factored out 
in one-photon absorption processes and approximately in the multiphoton case. This 
property allows to obtain cw perturbative photoionization cross sections from amplitudes 
extracted over the range of continuum energies within the spectral bandwidth of short 
laser pulses. We also study the formation of the profile of the Fano resonant peak 
in the time domain using ultrashort fs pulses, showing transient oscillations due to the 
two-path interference, which eventually vanish to yield the asymptotic stationary result. 

Atomic units are used unless otherwise stated. 



JPB 



3 



2. Theory 

The Feshbach projection method p2J [18] has shown to be a powerful method to 
describe resonance phenomena in scattering processes. Its application to the atomic 
electronic structure can be found elsewhere ([19] and references therein) although its 
practical implementation is mostly reduced to atomic systems with two and three 
electrons. A detailed study of the application of the stationary Feshbach method in 
He has been performed by Sanchez et al (20]. Also, after the pioneering work of 
Temkin and Bathia on three-electron systems [21], the Feshbach formalism has been 
recently revisited and applied to the Li atom in our laboratory [22J, including and 
assessing all the required ingredients of the rigorous formalism. The Feshbach projection 
operator formalism is based on the introduction of projection operators V and Q, 
satisfying completeness (V+Q=l), idempotency (V 2 =V, Q 2 =Q) and orthogonality 
(VQ=QV=0), which project the total wavefunction onto nonresonant scattering-like 
and bound-like halfspaces, respectively. These projected wave functions must also 
satisfy the asymptotic boundary conditions lim r .^ 00 V^/=^/ and lmv.^oo Q^/=0, the 
latter indicating the confined nature of the localized part of the resonance. 

By introducing the splitting of the total wave function ty=Q^ + into the time 
independent Schrodinger equation H^=Efy, working equations for the bound-like and 
the non-resonant scattering-like parts arise as follows (see, for instance, [22J): 

(QHQ-£ n )Q$ n = (la) 
(VH'V - E)V^° = 0, (lb) 

where H' is the operator containing the atomic Hamiltonian plus an optical potential 
devoid of any resonant contribution from the state <2$ s with energy S s , i.e., H'=H + 

Kpf ' where 

Vlt s = rf VHQ ^^ QHV. (2) 



opt ■ — ~~ p _ c 

It is worth noting that the Hamiltonian splits into H= QH Q+VHV+ QHV+VHQ, 
where the last two terms are responsible for the coupling between both halfspaces which 
ultimately causes the resonant decay into the continuum. In practice one starts by 
solving Eq. ( Ilaj) for the Q space with a configuration interaction method to obtain 
a first approximation to the location of resonant states and it implies to use Q—l-V, 
where V= V\+V2-V{P2 with Vi being a one-particle projection operator. In this work 
we are restricted to doubly excited states lying below the second ionization threshold 
of the He atom, so that Vi = |</>i s (rj))((/> ls (rj)|. Therefore, the Q operator removes all 
those configurations containing the Is orbital, then avoiding the variational collapse 
to the ground state (Is 2 ), to singly excited states (lsn£) and removing also the single 
ionization continuum (lsei). Instead of diagonalizing the QHQ problem, one may 
solve the equivalent but simpler eigenvalue problem involving an effective Hamiltonian 
that consist of the full Hamiltonian and a Phillips-Kleinman pseudopotential [23J, 
H e ff=H+MV, where M is a large real number. The effect of MV is to project 
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upward in energy all eigenstates associated with the V halfspace, so that the lowest 
variational eigenstates correspond to the Rydberg series of resonant discrete states below 
the corresponding ionization threshold. 

To build the nonresonant electronic continuum we solve Eq. (1 lb j) for the V 
subspace using a basis of two-electron configurations {^k}k=i, with cufc(xi,x 2 ) = 
A (</>i s (xx)^^)), where we use the notation Xfc=(rfc, Q^, s^) (radial r^, angular Q^, 
and spin Sk coordinates). This conforms a static- exchange approximation for the 
nonresonant continuum wave function V^°, which is orthogonal to Q, as proved by 
construction, since 0i s is used to build the V projection operator. As a general 
rule, the ground state, the singly and doubly excited states are obtained with a 
configuration interaction method, using a truncated set of antisymmetrized two- 
electron configurations {w 1 ^ L ' s }n=i adapted to the L, S total angular momenta 
and generated with products of atomic orbitals <p n e, as follows 

w l,m l ,s,m s{xi ^ 2) = A (0 ni , 1 (rO^ 2 (r 2 )^(^i^ 2 )x 5 ' Ms (^i,^)) , (3) 

where y^'^ 1 is the bipolar spherical harmonics, x the spin wave function, and A 
corresponds to the antisymmetrizing operator. In the stationary implementation of 
the Feshbach formalism, one may readily compute the energy shift correction A s for the 
resonant state Q$ s due to the perturbation of the surrounding continuum with 

As = ^ dE ,mQHPW))\\ (4) 

7^YE S E s - E' 

so that its actual energy is E s = £ s + A s . QHV couplings also allow to compute the 
resonant widths using a Fermi's golden rule, 

T s = 2tt\(^ s \QHP\V^°(E = E s )\ 2 (5) 



as well as the g-Fano profile parameter [20]. We have described above the basic 
ingredients to perform stationary calculations within the Feshbach formalism. We now 
proceed to describe our implementation of the time-dependent Feshbach method to be 
applied to atomic photoionization phenomena with ultrashort laser pulses. 

2.1. Time dependent Feshbach method. 

The time dependent Schrodinger equation (TDSE) for a two-electron atomic system 
exposed to a laser field in the dipolar approximation reads 

# + ^)-zJj*( Xl ,x 2 ,*) = (6) 

where ^(t) = (pi+p2)-A(i) (V r (t) = (ri+r 2 )-E(t)) corresponds to the laser-atom interac- 
tion in the velocity (length) gauge. The vector potential for a z-linearly polarized 
laser pulse is taken as A(t)=e z f(t)cos(oj(t — T/2)), where f(t) is a shape function for 
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the pulse envelope taken here as /(t)=sin 2 (7rt/T) and it is defined in the time inter- 
val t G [0,T] and zero elsewhere. The field amplitude is related to the laser intensity 
through /[W/cm 2 ] = 3.5095x 10 16 -Eq [a.u.] and to the vector potential amplitude with 
A = E /u. 

We make use of a spectral method by expanding the time-dependent wavefunction 
with the stationary Feshbach eigenstates of the QHQ and VHV projected Hamiltonians. 
This means that our expansion is based on the asymptotic Hamiltonian, since QHV 
couplings vanish not only for ri — > oo, but they are also not effective for t — > oo since 
eventually all the Q resonant population leaks to the continuum V halfspace. To work 
with an asymptotic basis of eigenstates has the advantage of avoiding projections of a 
propagated full wavepacket (without Feshbach partitioning) and therefore the calculated 
expansion coefficients Ci(t) for t ^> {T,l/T s } correspond to physical amplitudes by 
themselves. Then our expansion reads 

*(x ls x 2 ) = ^C b (t)S b ( Xl ,x 2 )e-^ 

b 

+ ^C r (t)Q$r(x 1 ,x 2 )e-^' 

r 

+ J dEC E (t)V^ E (*i^2)e- lEt (7) 

where {E b , Ef,}^ correspond to the set of bound states (ground plus singly excited 
states), {Q§ r , E r }^=i to the resonant doubly excited states and {V^° E , E}e refer to 
the nonresonant continuum states. Although rigorously the bound states pertain to 
the V space, in our implementation they are calculated separately as eigenstates of the 
full Hamiltonian H. Indeed, the CI basis of the static exchange approximation used 
to expand the V halfspace provides rather poor results for the bound states. Using 
E& states introduces non-orthogonalities in the basis set of configurations among states 
with the same L, S, n symmetries. For one-photon absorption (L — > L ± 1, S, ir — > it'), 
the lack of orthogonality is irrelevant, but for multiphoton processes it must be checked 
out carefully, introducing the corresponding overlaps in the couplings described below. 
Now introducing the ansatz of Eq. (J7|) in the TDSE (J6]) one arrives to a set of coupled 
differential equations in the interaction picture, which may be written in packed block 
form, with n={b,r, E}, as follows 
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where V{t) n%m are coupling matrix elements corresponding to the laser interaction and 
QHV are electrostatic couplings since neither Q$ nor are eigenstates of the 

field free Hamiltonian. The latter couplings do not depend on time and therefore, 
once the laser pulse ends, they are still active to empty the resonant population 
until its full depletion. The set of coupled equations is solved subject to an initial 
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condition C n (t=0)=5 n0 . In conclusion, the time-dependent Feshbach method relies on 
the preliminar calculation of the stationary eigenstates (QHQ and VHV) and couplings 
QHV of the projected Hamiltonian, along with the dipolar couplings among them. 

In other works (see, for instance, [2U [251 [26j [27J ) the total wavefunction is expanded 
in terms of the eigenfunctions of the complex-rotated (non-hermitian) unperturbed 
Hamiltonian H(re l6 ). Due to the presence of a negative imaginary part in the 
eigenenergies of both the resonant and the continuum states, the time propagation 
allows for the decay of the metastable state as well as to avoid non-physical reflections 
of the wave packet at the boundary box. Also, in the complex scaling (CS) method 
the complex eigenstate corresponding to the resonant pole includes both the bound-like 
and the scattering-like part of the resonant state and Fano profiles in photoionization 
spectra may be obtained straightforwardly with a reduced number of final complex 
eigenstates since each single continuum pole represents indeed a bunch of unrotated 
continuum states, i.e., takes into account the presence of neighbouring continuum 
levels. The description in terms of the Feshbach formalism requires instead a higher 
density of continuum states. The Feshbach approach treats separately the bound- 
like and the scattering-like parts of resonances, allowing for the population transfer 
between both halfspaces through the specific inclusion of the electrostatic couplings 
QHV. Concerning the laser-atom interaction, the dipolar operator in the CS method 
involves the simultaneous coupling of both the bound and the scattering part of the 
resonance with other states, let them be bound, continuum or another resonant state. 
At variance, in the Feshbach method the dipolar couplings are performed separately 
between wave functions located in V and/or Q halfspaces. 

Finally, multiphoton ionization cross sections can be obtained in the time domain 
using weak laser pulses with finite (but long enough) time duration (to be compared 
with time-independent perturbative results in the energy domain) using the expression 
(see, for instance, [T5] ) 

r 27V TV-li _ ( u[Joules} \ " P TDS E {) 

a[Cm S l ~ \l[W/cw?\) C{N)[s\ W 

where u is the central frequency of the absorbed laser pulse, / the laser intensity and 
C(N)=Jq dtf(t) 2N takes into account the effective pulse duration using its envelope 
f(t) and it consists of a numerical factor times the pulse duration, i.e., C(1)=3/8T, 
C(2)=35/128T, and so on. The total ionization probability is directly obtained 
from the coefficients Ce in the numerical solution of the coupled equations ([H]) with 
PTDSE=Yl l E i \^Ei(t > T)\ 2 , i.e., just a summation over all discretized continuum states 
(a simple quadrature of the integral) [28J. 
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3. Computational details 

3.1. Bound states. 

Our CI method is based on the expansion in terms of antisymmetrized products of 
atomic orbitals, the latter expanded in B-splines polynomials enclosed within a finite 
box of length L. B-splines have been widely used in the last years and for a fuller 
description the reader is referred to [29J. A very precise ground state energy for He 
atom can be obtained with B-splines basis using an exponential knot sequence and 40 
B-splines with order k—7, generating one-electron orbitals with I < 3, with a full CI wave 
function of 3280 configurations, which yields the energy -2.903321 a.u. to be compared 
with that of Pekeris [30], -2.903742 a.u. Nevertheless, since one is mainly interested 
in processes where the electronic continuum is mostly involved, the exponential knot 
sequence is not well adapted to the description of continuum wave functions, due to the 
poor electronic density in the energy region of interest. Instead, we have used in this 
work an exponential- linear sequence of knot-points within a box of length L=150 a.u., 
with 200 B-splines, in such a way that the sequence is exponential below 15 a.u. and 
linear in the rest of the box. This allows for a good description of the oscillations of 
the continuum wave functions in the full range of the box. Therefore, 27 bound states 
of symmetry 1 S e (10 states), 1 P° (9) and 1 D e (8) are obtained with 8230, 9248, and 
10986 configurations, respectively, using two-electron angular configurations built with 
orbitals with £ < 3. Eventually, for the 1 S e ground state, the basis of orbitals may be 
supplemented with Slater-type orbitals expanded in terms of B-splines to improve the 
convergence of the ground state energy [201 EI] • 

3.2. Q-subspace. 

For the QHQ resonant space, we perform CI calculations with the same configurational 
basis set but removing the Is orbital. Then we are able to obtain 19 doubly excited 
states of symmetry 1 S e , 26 states for 1 P°, and 25 states for 1 D e , using 8456, 9135 and 
10861 configurations, respectively. We do not intend here to produce benchmark results, 
since we have not followed a systematic optimization procedure, but to produce fairly 
good results to show the performance of the time- dependent Feshbach method. To 
illustrate the accuracy of our computation, we report in table CD the energies (including 
the energy shift correction of Eq. (j3j)) and widths (or lifetimes) for the 1 S e , 1 P° and 
1 D e lowest doubly excited states below the N=2 threshold, and compared to accurate 
results obtained by Chen [32] using the saddle-point complex rotation method. For a 
more comprehensive comparison, the reader is referred to [32]. 

3.3. V-subspace. 

Since the illustration of the method is restricted to one- and two-photon ionization 
below the N=2 threshold in He atom, we only deal with the simple case of single- 
channel continua. The multichannel case for processes above N=2 within the 
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Table 1. Energy positions (in a.u.), witdhs (in a.u.) and lifetimes (in fs unless 
otherwise indicated) for the lowest four 1 S e , 1 P° and 1 D e doubly excited states in He 
1 S e below the N=2 threshold. The values of the present work are compared with the 
saddle-point complex- rotation results of Ref. [32]. The notation a[b] indicates axlO b . 



Sym. 


State 


£ r +A s 


r s 


r [fs] 




2(1,0) + 


-0.777533 


0.5051[-2] 


4.79 




m 


-0.77787 


0.453[-2] 






2 (-i.o) + 


-0.619822 


0.2419[-3] 


100.0 




[32] 


-0.62181 


0.2178[-3] 






2 (i, o)t 


-0.589631 


0.1507[-2] 


16.1 




[32] 


-0.589896 


0.137[-2] 






2 (-i,o) + 


-0.547755 


0.9895[-4] 


244 




[32] 


-0.548070 


0.775 [-4] 




1 po 


2(0, l)t 


-0.692642 


0.1392[-2] 


17.4 




[32] 


-0.693069 


0.1372[-2] 






2(1,0)3 


-0.597065 


0.4069[-5] 


5945 




[32] 


-0.597074 


0.384[-5] 






2(0, i)t 


-0.563721 


0.2956[-3] 


81.8 




[32] 


-0.564074 


0.2998 [-3] 






a(-l,0)g 


-0.547031 


0.2358[-7] 


> 1 ns 




[32] 


-0.547087 


0.15[-7] 




l D e 


2(1,0)+ 


-0.701512 


0.2528[-2] 


9.56 




[32] 


-0.70183 


0.236[-2] 






2(1, o)i 


-0.569448 


0.6052[-3] 


39.9 




[32] 


-0.569193 


0.560[-3] 






a(0,l)S 


-0.556317 


0.2015[-4] 


1200 




[32] 


-0.556417 


0.201 [-4] 






2(1,0) 4 + 


-0.536575 


0.2485[-3] 


97.3 




[32] 


-0.536715 


0.234[-3] 





Feshbach formalism has also been developed [201 I2H1 [3TJ [33]. A set of discretized 
nonresonant continuum states with energies E\ and angular momentum L are 

obtained by diagonalizing Eq. ( ffbl for the V space using a CI expansion V^l\ = 
J2f=i CjA((fii s ■ (pj£=L) where one of the electrons is frozen in the </> ls orbital (static 
exchange approximation) and the other <ftj£ pertains to the subset of iV hydrogenic 
orbitals generated with M B-splines with angular momentum £=L. Therefore, the 
wave function for the ejected electron corresponding to the state with discretized energy 
£'i>-2.0 (first ionization threshold) can be built using the CI expansion coefficients 
[31], i-e., ipi(r)— $2 3 -=i Cj^W- This discretized continuum wavefunction is normalized 
to unity since it comes from a diagonalization procedure. In order to renormalize it 
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to the correct Dirac delta, it suffices to multiply with the density of states factor, 
i.e., -ipE i = [p(Ei)} 1 ^ 2 'ipi, computed with a two-point formula from the set of discretized 
continuum energies, p(Ei) = 2/(E i+ i — [35J. The function r • ip Ei contains 

the effect of the Coulomb and the short-range potentials, and its radial part behaves 
asymptotically as ^2/Trk Ei \Fki{k Ei ,r) cos 5 El + Gki(k Ei ,r) sin 5^] [36], where Tkt and 
Qki correspond to the regular and irregular Coulomb functions, respectively, and 5 Ei 
is the scattering phase shift against the free Coulomb wave. By a least-square fitting 
procedure the computed function r ■ ip E . can be adjusted to the analytical asymptotic 
form inside an interval [r asym , L] in the outer part of the box of size L (covering at least 
two wavelengths), from which the corresponding phase shifts 5 Ei can be obtained. 

The knowledge of these phase shifts is useful not only to compute differential 
cross sections, but also to provide a prescription to adjust the discretized continuum 
energies {E{\f =l to any desired energy value. For instance, to compute the widths with 
Eq. dSJ), the nonresonant continuum state must be degenerate with the resonant 
energy E s . By diagonalization of the V eigenvalue problem, none of the Ei energies 
coincide with E s , but E s interpolates two neighboring discretized continuum energies, 
Ei < E s < E i+ \ and so do the phase shifts S Ei < 5 Es < S Ei+1 . Then, the interpolated 5 Es 
is introduced in the asymptotic radial formula given above. Every single-electron radial 
continuum function obtained by diagonalization inside the box satisfy the boundary 
condition r • tp Ei (r=L)=0, but the asymptotic analytical form with the interpolated 
phase shift 5 E does not. Instead, it shows a node at r very close to the edge of the box 
L. This means that in order to get one of the eigenvalues Ei be equal to E=E S the box 
length should be replaced by r . To avoid changes in the box length and therefore in 
the B-splines basis itself, one may keep the same box length L, but with the addition of 
a step potential V Q Q{r — r ) (0 is the Heaviside function and V is a large real number) 
in the eigenproblem for the atomic orbitals. This prescription guarantees that at least 
one of the eigenvalues Ei is nearly degenerate with the selected E value. This general 
procedure can be applied to the computation of resonance widths as well as to produce 
a new set of discretized continuum states {ei}fii N with an energy spacing selected at 
will, which may be useful, for example, to increase the density of states if required. 

In this work we have computed the nonresonant continuum states for 1 S e , l P° and 
1 D e (L=0,l,2) with CI ((j)u<f>ji=L) configurations, where j runs up to 170 orbitals for 
each £=L angular momentum. At least, 75 continuum states are obtained lying above 
the first ionization threshold and below the second one, the energy region where the 
lowest doubly excited states of Q-space are located. 

S.Jf.. Time- dependent calculations 

Once the static calculations for the Q and V space are performed, the energies, dipolar 
couplings and the QHV electrostatic couplings must be introduced in Eq. (jSJ), a 
system of differential coupled equations which is solved using a sixth-order Runge-Kutta 
integrator. Long time propagations (even beyond the laser pulse duration) require large 
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Figure 1. (color online) Time decay of the three lowest 1 P° resonances in He, 2(0, 1)J 
(blue), 2(1,0)3 (red) and 2(0, 1)^ (green) quoted in table [TJ obtained by solving the 
field-free TDSE with only QHV couplings, for a propagation time of t=100 fs, and 
compared to the exponential decay formula P(t)=P(0)e~ r ' * (dotted lines), where T r 
values are taken from tabic [1] 

radial boxes and energy spacings in the discretized nonresonant continuum smaller than 
the spectral width of the pulse (approximately given by Au= Att/T) and at least smaller 
than the largest resonance widths T s . With only 75 continuum states between the first 
and the second ionization threshold, the density of states is not large enough, but one 
may create a larger boxby simply interpolating the dipolar and QHV couplings involving 
the continuum with energies {Ei]f =l . A new much denser set of continuum energies 
{E'j}^ can be generated following a quadratic formula E'—Eo+ai 2 , with a=(e%— e )/M 2 
where Eq and E\ correspond to the energies of the first and second ionization thresholds, 
and M is the number of interpolating grid points. In our calculations we have used 
up to 2000 interpolating points for the dipolar and QHV couplings. This method 
does not fully replace the consistency of using larger boxes but for the purposes of our 
calculations it serves in practice to avoid unphysical spreading of the wavepacket and 
undesired reflections from the box boundary, for the propagation times used in this work. 
For the results reported in the next section we always include all the computed bound 
and resonant states, and the set of interpolated continuum states for three symmetries, 
1 S e , l P° and l D e , a minimal angular basis for two-photon transitions. 

4. Results. 

To illustrate the performance of the time-dependent Feshbach method in He, we 
restrict ourselves to the study of one- and two-photon ionization with laser intensities 
corresponding to the perturbative regime, and moderately long femtosecond pulses. 
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Time [fs] 

Figure 2. (color online) Time resolved photoexcitation and decay of the two fastest 
decaying 1 P° resonances in He, 2(0, 1)2 (upper panel) and 2(0, 1)3" (lower panel) in the 
one-photon ionization process from the ground state, obtained by solving the TDSE 
for He subject to a laser field with 7=10 10 W/cm 2 and duration T—20 fs, with a total 
propagation time t= 100 fs. The color scheme corresponds to computations with central 
frequencies u corresponding to different small positive (red) and negative (black) 
detunings, the line in green with the largest population corresponds to zero detuning. 
The dotted lines indicate the exponential decay according to P(t)=P(T)e~ rl - t ~ T \ 

4.1. One-photon ionization from the 1 S e ground state 

First, in order to test the adequacy and effectiveness of the computed QHV couplings, 
responsible for the decay of the doubly excited states into the continuum, we consider 
field-free time propagations of the TDSE up to 100 fs, in which the initial state 
corresponds to a given fully populated resonance, i.e., C r (£=0) = l. This should represent 
a non stationary evolution, and the expected time-dependent probability must follow 
an exponential decay P r (t)=P r (0)e~ Trt , where T r corresponds to the resonance width, 
computed with the Feshbach stationary method (see figured]). 

According to table [1] only the first and the third lowest 1 P° doubly excited 
states have lifetimes under 100 fs, and therefore within this limit of time propagation 
only features corresponding to these two resonances would be noticeable in the 
photoionization spectra. One-photon ionization cross section from the ground state 
1 S e + wH' 1 ?" is computed by solving the TDSE with both dipolar and QHV couplings 
included, for a laser pulse with intensity i=10 10 W/cm 2 , duration T=20 fs and 
propagation time up to £=100 fs, and using a grid of photon energies 24.5 eV < E u < 
68 eV. For each photon energy, the total ionization probability PrDSEif) is collected 
at times when the amplitude of the vector potential A(t) vanishes (the nodes of the 
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Figure 3. (color online) Temporal built-up of the one-photon ionization cross section 
from the ground state of He to the 1 P° continuum between the first He + (7V=1) and 
second He + (N—2) ionization threshold. Cross sections according to Eq. ^ are given 
in Mbarn, photon energies in eV and propagation time in femtoseconds. Length and 
velocity gauge results are indistinguishable in the figure. 



function chosen for A(i)), to comply with a proper field- free projection during the 
laser pulse. Ionization probabilities during and after the pulse duration come from 
interfering amplitudes due to the laser direct ionization and to the time delayed non 
stationary decay of resonances, and this interference is known to produce the typical 
Fano profiles in the photoabsorption spectra. When the laser field is present, and 
for photon energies close to the resonant condition E u = E r — ^(l 1 ^), the doubly 
excited states are populated from the ground state from the onset of the laser pulse, 
reaching the maximum population in the second half of the pulse duration, from which 
the resonances manifestly decay into the continuum with the exponential law. The 
time-resolved photoabsorption and time-delayed decay of the two most representative 
(fast-decaying) 1 P° resonances is pictured in figure |2] for different photon detunings. 

Fano interferences take place during the laser pulse but only at very long 
propagation times (longer than the lifetime of the slowest decaying resonance in 
the Rydberg series) the Fano profile reaches its asymptotic form for all resonances. 
Nevertheless, the onset and time evolution of the Fano profile corresponding to the 
fastest decaying 1 P° resonances can be studied within the practical limitations of the 
present application due to the limited chosen basis and box size. In figure [3] we plot 
the time dependent built-up of the photoionization cross section between the first and 
the second ionization thresholds for the same laser pulse described above. In fact, a 
laser pulse with duration T=20 fs allows to visualize the formation of the nonresonant 
background as well as the two major resonant peaks present in the time-independent 
perturbative spectrum. After the laser pulse duration there are slight changes in the 
total ionization and after t=50 fs the behavior is almost stationary. From this figure we 
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Figure 4. (color online) Time-resolved formation of the Fano profile in the differential 
photoionization probability dP/dE in the energy region of the lowest He 1 P° resonant 
state, using a laser pulse with central frequency w=2.21 a.u., intensity 7=10 10 W/cm 2 
and duration T=20 fs (the time evolution of the laser pulse is also drawn). Snapshots 
taken at 5, 10, 20, 30 and 50 fs. The red solid line corresponds to the result obtained 
with an analytical model [7] (see text). The blue dots are calculated with the time- 
dependent Feshbach formalism. 



may appreciate that the resonant peaks appear very soon after the pulse is switched on, 
but their profile is almost symmetrical. Only in the second half of the laser pulse, the 
resonant states become noticeably populated, subsequently decaying and the interfering 
with the direct ionization, which produces ultimately the final asymmetry in the profiles. 

The time dependent formation of the Fano profile can also be analyzed from the 
photoionization probabilities, differential in energy, dP(E,t)/dE. To obtain reliable 
differential probabilities in the continuum a rather good density of continuum states 
close to the resonance positions is required, and they are simply computed with the 
continuum expansion coefficients, i.e., dP(E,t)/dE = p(E)\C Eif)\ 2 , where p(E) is the 
density of states. Nicolaides et al [7] have also recently studied the time dependent 
formation of the profile of the lowest He 2s2p (or 2 (0, 1)^) 1 P° resonant state, induced 
by a short laser pulse. Apart from their ab initio calculation using an spectral method 
called state specific expansion approach, they propose a simplified theoretical model to 
account for the profile formation. In their simplified model the transition amplitude to 
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Figure 5. (color online) The same as in figure |4] but for the third lowest He 1 P° 
resonant state, then using a laser pulse with central frequency w=2.34 a.u. 



the continuum is given by 

A E (t) = C E e~ iEt G(E,t) 

- C E ^e-^ +Ar - ir ^G(£ r + A r - iT r ,t), (10) 

where {£ r , A r , T r , q} are the set of resonance parameters corresponding to uncorrected 
position, energy shift, width and q Fano shape parameter, respectively, the coefficient 
C e=F r d g E / — (£~ r + T r )] 2 + T^, with d g E being the dipolar coupling value between 
the ground state and the scattering state with energy E. The G(E, t) function contains 
the dependence on the laser pulse and it is basically the Fourier transform of the electric 
field of the laser pulse, i.e., G(E, t)=—i J * dt' e~ li * Eg ~ ES>t ' ~E,{t)e z . We compare our ab initio 
results with this simplified model, for the same laser parameters described above with a 
central frequency w=2.211 a.u., J=10 10 W/cm 2 and duration T=20 fs and the resonance 
parameters quoted in table [[J with q=-2.8 and d g E=0. 48185. We have computed ab 
initio differential ionization probabilities for continuum energies close to the lowest 1 P° 
resonance, -0.75 < _E[a.u.] < -0.63, and they are plotted in figure HI We notice that the 
analytical model reaches the asymmetric profile much more rapidly than the ab initio 
method, the latter showing the effect of laser spectral width at short times, in the sense 
that the He atom is not yet aware of the total pulse duration. The Feshbach results only 
converge to the expected asymptotic result of the model for times much longer than the 
pulse duration, with small transient oscillations due to interferences associated to the 
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Figure 6. (color online) Top panel: One-photon single ionization cross section 
[in Megabarns] , against the photon energy [in eV] of He, obtained solving the 
Feshbach TDSE with a laser pulse of duration T=20 fs and intensity 7=10 12 
W/cm 2 , with a propagation time of £=100 fs. Middle panel: One-photon single 
ionization cross sections, extracted from the differential photoionization probabilities 
dP/dE=p(E)\CE\ 2 plotted in the bottom panel, obtained using a laser pulse with 
shorter duration T=0.9 fs and intensity 7=10 12 W/cm 2 and central frequencies 
w=27.57 eV (black), 36.76 eV (red), 45.95 eV (green), 55.14 eV (blue). Sharper 
signals for the resonant peaks are obtained extracting cross sections from differential 
probabilities (scaled with a factor 10 -2 in the figure) calculated with longer pulse 
duration T=10 fs for photon energies between ^57 and >~65 eV. Vertical dashed 
lines indicate ionization thresholds He + (N=1) and He + (iV=2). All reported results in 
velocity gauge. 

still active resonance decay. 

The same comparison is done with the third lowest 1 P° resonance in figure [5] in the 
energy range -0.62 < J5[a.u.] < -0.5, now with a central frequency u=2.34 a.u. a dipolar 
coupling value d g E=0. 439891 a.u. and a Fano parameter q=-2.5, and the resonance 
parameters quoted in table [TJ Again, the Feshbach result nicely shows the trend of 
convergence to the analytical result of the model for propagation times much larger 
than the pulse duration. Of course, our ab initio simulations in the present illustrations 
are constrained by our box length L=150 a.u., and therefore the total propagation 
time is also limited. In fact, this drawback in time dependent computations has been 
recently overcome. Usually, autoionizing lifetimes may be much longer than the pulse 
duration and the required field-free integration times are prohibitive in order to obtain 
asymptotic ionization probabilities in the neighborhood of all resonances conforming 
a Rydberg series. In a recent series of papers by Palacios et al [37J EH] it is shown 
that multiphoton ionization cross sections can be retrieved from differential amplitudes 
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obtained with very short laser pulses, the asymptotic amplitudes being extracted by 
solving a driven equation within the exterior complex scaling approach. This method is 
in practice equivalent to propagate the time-dependent wave-packet for an infinite time 
after the pulse duration. This methodology could be eventually incorporated to our 
Feshbach time-dependent method, but this will not be done here. Instead, we use their 
procedure to obtain photoionization cross sections (for which one assumes low intensities 
and laser pulses with infinite duration) from differential ionization amplitudes obtained 
with short laser pulses, thanks to the factorability of the time dependence in first-order 
perturbation theory, in the form of the Fourier transform of the pulse. 

For instance, using expression (17) in Ref. [38], modified to be used with our 
differential probability amplitudes in the energy scale, we can reproduce both the 
nonresonant background and the resonant sharp peaks in one-photon ionization cross 
section using transition probabilites obtained within the spectral bandwidth of shorter 
laser pulses with durations T=0.9 fs and T—W fs, the latter to improve the resolution 
in the resonance region (see figure [6]). Again, since we are not using truly asymptotic 
amplitudes (t — > oo) we cannot obtain the high resolution in Fig. 2 of Ref. [38J for the 
Rydberg series of 1 P° resonances. 

4-2. One-photon ionization from the lowest ls2p 1 P° excited state. 

In order to reach also the lowest 1 S e and 1 D e resonances, one may check out the one- 
photon ionization from the ls2p 1 P° state. According to tabled] the lowest resonances 
in these two symmetries decay even faster than those in l P°. The cross section obtained 
for a laser pulse of duration T=20 fs and total propagation £=100 fs is plotted in figure CD 
The observed peaks in the spectrum correspond to the three lowest 1 S e doubly excited 
states quoted in table [1] and the lowest 1 D e resonant state. Our cross sections compare 
well with the perturbative stationary result, computed with a Multiconfigurational 
Hartree-Fock (MCHF) method and B-splines [39J. The ls2p 1 P° excited state is closer 
to the upper ionization than the ground 1 S e state. A photon energy of ~21.2 eV is able 
to photoionize the He atom, but this energy is also resonant with the ground 1 S e state. 
In this case the laser field can couple different one- and three-photon processes to yield 
1 S e and l D e enhanced ionization probabilities, along with Rabi oscillations between 
the initial state l P° and the ground state 1 S e , nonlinear phenomena hardly seen at this 
perturbative low laser intensities, but expected to play a major role at higher intensities. 

4-3. Two-photon ionization from the 1 S e ground state. 

For the sake of completeness of this dynamical study, we also compute the two-photon 
ionization process from the ground state, in which the 1 S e and 1 D e doubly excited 
states can be also excited with the absorption of a first photon into the intermediate 
1 P° continuum (above threshold ionization). Here, continuum-continuum dipolar matrix 
elements are explicitly included, and in this case it is known that calculations for above 
threshold ionization in the velocity gauge converge faster that in length gauge [40J, 
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Figure 7. (color online) One-photon ionization cross section from the \s2p 1 P° 
excited state to partial final contributions 1 S e (red) and 1 D e (green), obtained from our 
Feshbach time-dependent method, using a laser pulse with intensity /=10 10 W/cm 2 
and duration T=20 fs. Peaks denoted with S correspond to 1 S e doubly excited states 
and that with D correspond to the lowest 1 D e doubly excited state, all quoted in Table 
[T] Results in velocity gauge are only included. 

so our results in figure |8] are given only in the velocity gauge. A first computation 
(top panel in figure E]) is carried out with a laser pulse with intensity 7=10 10 W/cm 2 
and a pulse duration T=20 fs, with a total integration time t=100 fs and we plot 
the separate x S e and l D e partial contributions. The peaks corresponding to resonant 
intermediate 1 P° bound states in the photon energy region [21.2,24.6] eV are reasonably 
well resolved. In the above threshold ionization region u > 24.6 eV, the two major 
features correspond to peaks generated by the 1 S e (2s 2 or 2 (1, 0) J) lowest resonant state 
and the l D e (2p 2 or 2 (1, 0)^) lowest state. A better resolution for the resonant peaks in 
the above threshold ionization region can be achieved by extracting the cross sections 
from the photoionization amplitudes obtained by a short pulse and renormalized by the 
two-photon shape function of the laser pulse (see Ref. [38] for details). 

Adapting the expression (25) in [38] to the energy scale and using our 
photoionization amplitudes p 1 ^ 2 (E)Ce (which in our case are not truly asymptotic) we 
may reproduce the two-photon ionization cross sections piecewise, using an assortment 
of laser pulses of duration T=5 fs with central frequencies from 12.4 to 32.26 eV. In 
panel A of figure [8] we show the piecewise reconstruction of the 1 S e component in the 
total cross section, obtained from differential ionization probabilities p(E)\Ce\ 2 shown 
in the panel B below. Similarly, the 1 D e component in panel C is generated following 
the same procedure with the differential probabilities shown in panel D. With this 
pulse renormalization procedure (valid only in the perturbative regime), the structure 
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Figure 8. (color online) Two-photon ionization cross section from the 1 S e ground 
state in He. Top panel: 1 S e (red) and 1 D e (green) components of the cross section, 
computed by solving the Feshbach-TDSE with pulses of intensity 7=10 10 W/cm 2 
and duration T—20 is. Panel A: 1 S e component of the two-photon cross section 
applying the piecewise renormalization procedure of Ref. [38] and using the computed 
differential probabilities included in panel B, obtained with laser pulses with duration 
T=5 fs and the same intensity. Every segment in the reconstructed cross section is 
shown with a different color scheme. Panel C an D: The same as Panel A and B, but 
for the 1 D e component of the total cross section. The dashed vertical line indicates 
the above threshold ionization threshold. 



of resonant peaks in the Rydberg series are more cleanly discriminated, but at the cost 
of partially loosing the structures due to one-photon absorption to intermediate 1 P° 
bound states, which require longer laser pulse durations to be consistently resolved [38J . 



5. Conclusion 



In conclusion, this work describes the theoretical details and the inner workings of an 
ah initio time dependent Feshbach method as applied to the resonant photoionization of 
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He atom using ultrashort laser pulses, below the second ionization threshold. Some 
simple illustrations of the performance of the method have been included, related 
to one and two-photon ionization processes, without further sophistications. It is 
assumed that there is much room for improvements, in terms of a optimized grids 
of B-splines, much larger radial boxes, larger configurational basis and partial waves for 
better convergence, the introduction of complex absorbing potentials, or even exploring 
other time propagators. At the present level, all computations presented here have 
been performed in simple desktop computers. In order to deal with pump-probe-like 
experiments where a XUV pump laser excites the resonances to be subsequently probed 
by an IR laser field, the multichannel extension of the method is required, but it can be 
implemented straightforwardly within the Feshbach formalism [29J. Steps along these 
directions are under way. 
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